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Abstract 

Populations of species in ecosystems are often constrained by availability of resources within their 
environment. In effect this means that a growth of one population, needs to be balanced by com¬ 
parable reduction in populations of others. In neutral models of biodiversity all populations are 
assumed to change incrementally due to stochastic births and deaths of individuals. Here we pro¬ 
pose and model another redistribution mechanism driven by abrupt and severe collapses of the entire 
population of a single species freeing up resources for the remaining ones. This mechanism may be 
relevant e.g. for communities of bacteria, with strain-specific collapses caused e.g. by invading 
bacteriophages, or for other ecosystems where infectious diseases play an important role. 

The emergent dynamics of our system is cyclic “diversity waves” triggered by collapses of globally 
dominating populations. The population diversity peaks at the beginning of each wave and expo¬ 
nentially decreases afterwards. Species abundances are characterized by a bimodal time-aggregated 
distribution with the lower peak formed by populations of recently collapsed or newly introduced 
species, while the upper peak - species that has not yet collapsed in the current wave. In most 
waves both upper and lower peaks are composed of several smaller peaks. This self-organized hier¬ 
archical peak structure has a long-term memory transmitted across several waves. It gives rise to 
a scale-free tail of the time-aggregated population distribution with a universal exponent of 1.7. 
We show that diversity wave dynamics is robust with respect to variations in the rules of our model 
such as diffusion between multiple environments, species-specific growth and extinction rates, and 
bet-hedging strategies. 


Author Summary 

The rate of unlimited exponential growth is tradition¬ 
ally used to quantify fitness of species or success of or¬ 
ganizations in biological and economic context respec¬ 
tively. However, even modest population growth quickly 
saturates any environment. Subsequent resource redis¬ 
tribution between the surviving populations is assumed 
to be driven by incremental changes due to stochastic 
births and deaths of individuals. Here we propose and 
model another redistribution mechanism driven by sud¬ 
den and severe collapses of entire populations freeing up 
resources for the growth of others. The emergent prop¬ 
erty of this type of dynamics are cyclic “diversity waves” 
each triggered by a collapse of the globally dominating 
population. Gradual extinctions of species within the 
current wave results in a scale-free time-aggregated dis¬ 
tribution of populations of most abundant species. Our 
study offers insights to population dynamics of microbial 
communities with local collapses caused e.g. by invading 
bacteriophages. It also provides a simplified dynamical 
description of market shares of companies competing in 


an economic sector with frequent rate of bankruptcy. 


Introduction 

Mathematical description of many processes in biology 
and economics or finance assumes long-term exponential 
growth [El yet no real-life environment allows growth 
to continue forever la. In biology any growing popula¬ 
tion eventually ends ups saturating the carrying capacity 
of its environment determined e.g. by nutrient availabil¬ 
ity. The same is true for economies where finite pool of 
new customers and/or natural resources inevitably sets 
a limit on growth of companies. Population dynamics in 
saturated environments is routinely described by neutral 
“community drift” models 6j sometimes with addi¬ 
tion of deterministic differences in efficiency of utilizing 
resources Q- 

Here we introduce and model the saturated-state dy¬ 
namics of populations exposed to episodic random col¬ 
lapses. All species are assumed to share the same envi¬ 
ronment that ultimately sets the limit to their exponen- 
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tial growth. Examples of such systems include local pop¬ 
ulations of either a single or multiple biological species 
competing for the same nutrient, companies competing to 
increase their market shares among a limited set of cus¬ 
tomers, etc. Specific case studies can be drawn from mi¬ 
crobial ecology where susceptible bacteria are routinely 
decimated by exposure to bacteriophages (see e.g. [s. 9] 
and references therein), or paleontological record, where 
entire taxonomic orders can be wiped out by sudden ex¬ 
tinctions happening at a rate independent of order size 

Q. 

Model 

Population growth P(t) of a single exponentially repli¬ 
cating species in a growth-limiting environment is tradi¬ 
tionally described by Verhulst’s [4| or logistic equation 
dP/dt = Q • P • (1 — P/Ptot), where the carrying capac¬ 
ity of the environment P tot without loss of generality can 
be set to 1. In this paper we consider the collective dy¬ 
namics of multiple populations competing for the same 
rate-limiting resource: 

• Local populations are subject to episodic random 
collapses or extinctions. The probability of an ex¬ 
tinction is assumed to be independent of the pop¬ 
ulation size. Immediately after each collapse the 
freed-up resources lead to the transient exponential 
population growth of all remaining populations Pi . 
The growth stops once the global population JT Pj 
reaches the carrying capacity P tot = 1. 

• The environment is periodically reseeded with new 
species starting from the same small population size 
7 < 1 (measured in units of environment’s carrying 
capacity). 

We initially assume that the growth rates and collapse 
probabilities of all species are equal to each other. We 
also disregard the neutral drift (Ej in sizes of individ¬ 
ual populations during the time between subsequent col¬ 
lapses. All these assumptions will be relaxed, simulated, 
and discussed later in the paper (see Supplementary Ma¬ 
terials SI Text, S1-S7 Figures). The number of species in 
the steady state of the model is determined by the com¬ 
petition between the constant rate of introduction of new 
species and the overall rate of extinctions in the environ¬ 
ment that is proportional to the number of species. To 
simplify our modeling we will consider a closely related 
variant of the model in which the number of species N is 
kept strictly constant. In this case each extinction event 
is immediately followed by the introduction of a brand 
new species. We have verified that the two version of 
our model have very similar steady state properties. The 
dynamics of the fixed-TV model is described by 

dPi/dt = Q-Pi-(l-J2 p j)-Vi(t)-Pi , ( 1 ) 
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where r\i (t) is the random variable which is equal to zero 
for surviving populations and has a large positive value 
for populations undergoing an extinction/collapse. 

To speed up our simulations we do not continuously 
calculate Eq. O since most of the time the carrying 
capacity of the environment is saturated when local pop¬ 
ulations do not change. Instead we simulate the model at 
discrete time steps marked by extinction events. At ev¬ 
ery time step a randomly selected local population goes 
extinct and a brand new species with population 7 <C 1 
is added to the environment. We then instantaneously 
bring the system to its the carrying capacity by multi¬ 
plying all populations by the same factor. 


Results 

In spite of its simplified description of the ecosystem 
disregarding pairwise interactions between species our 
model has a rich population dynamics. Figure [l]V shows 
the time-courses of populations in a system with N = 20 
species and 7 = 10 -4 . At certain times the carrying 
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FIG. 1. Population dynamics. The simulated model has 
N = 20 species and 7 = 10 -4 . (A) Time-courses of popula¬ 
tions of all species. The color denotes population sizes (see 
the color scale on the right) with the dominating species visi¬ 
ble as red horizontal bands. Note five diversity waves ending 
at purple dashed lines. Transitions between these waves were 
triggered by extinctions of the dominating species # 5, 15, 6 , 
19, 16 correspondingly. (B) The time-course of the species # 6 
with the logarithmic y-axis. Note the pattern of intermittent 
periods of exponential growth fueled by local extinctions. 

capacity of the environment is nearly exhausted by just 
one dominant species with P macc ~ 1 visible as dark red 
stripes in Figure [l]V. When such dominant species goes 
extinct a large fraction of the resources suddenly becomes 
available and consequently all other populations increase 
by a large ratio 1/(1 — Pmax )• This marks the end of 
one and the start of another diversity wave that initially 
is dominated by a large number of species with similar 
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population sizes. In the course of this new wave these 
species are eliminated one-by-one by random extinctions 
until only one dominant species is left standing. Its col¬ 
lapse ends the current and starts the new wave. In Figure 
UK one can clearly distinguish about 5 such waves termi¬ 
nated by the extinctions of dominant species #5, 15, 6 , 
19, and 16 correspondingly. 

Figure UK shows the time-course of just one local pop¬ 
ulation of the species #6 on a logarithmic scale. Be¬ 
tween time steps 100 and 150 the population of this 
species nearly exhausts the carrying capacity of the en¬ 
vironment. Its local extinction at the time step 154 
ended the third diversity wave and started the fourth one. 
Note somewhat erratic yet distinctly exponential growth 
of this species happening on the slow timescale set by 
the frequency of extinctions. This growth should not be 
confused with exponential re-population of recently col¬ 
lapsed environments that happens much faster (a small 
fraction of one time step). 

Figure [2] follows the population diversity (grey shaded 
area) defined as D(t) = 1/X^Li Pi(t) 2 as a function 
of time in a system of size N = 1000. In general D 
can vary between ~ 1 when one abundant species domi¬ 
nates the environment and N when all species are equally 
abundant. The diversity is inversely proportional to the 
largest population P m ax(t) = max^ Pi(t). The diversity 
waves (purple dashed arrows in Figure [ 2 ]) are initiated 
when a dominating species collapses. As a consequence, 
at the start of each wave the diversity abruptly increases 
from ~ 1 to a substantial fraction of the maximal pos¬ 
sible diversity N. After this initial burst triggered by 
the global redistribution of biomass, the diversity expo¬ 
nentially declines as exp (—t/N) (the dot-dashed line in 
Fig|2|), driven by ongoing extinctions of individual pop¬ 
ulations. The typical duration, t wave of a diversity wave 
is equal to the time required for all of the species present 
at the start of the wave to go extinct one-by-one. Thus 
it is determined by N • exp(— t waV e/N) ^ 1 or 

twave rv N • log e N . (2) 

Figure [3] shows the time-aggregated distribution of 
population sizes for 7 = 10 -9 and N = 1000 (the grey 
shaded area). This logarithmically-binned distribution 
defined by 7 r(P) = dProb(Pi(£) > P)/dlog 10 P were col¬ 
lected over 20 million individual collapses (time-steps 
in our model). Thus, a time-aggregated distribution is 
rather different from a typical “snapshot” of the system 
at a particular moment in time characterized by between 
1 and N of highly abundant species in the current diver¬ 
sity wave. The time-aggregated distribution is bimodal 
with clearly separable peaks corresponding to two popu¬ 
lation subgroups. The upper peak consists of the species 
that have not yet been eliminated in the current wave. 

To better understand the dynamics of the system in 
Figure [3] we also show the distribution of populations 
sizes at the very end of each diversity wave (green line) 
and at the beginning of the next wave (red line). That 
is to say, for the green curve we take a snapshot of all 



FIG. 2 . Diversity dynamics. The grey shaded area shows 
the the time course of the population diversity D — 1/ JT Pf 
in our model with N = 1000 and 7 = 10 -12 . Purple dashed 
lines mark the beginnings of diversity waves when a collapse 
of the dominant species with Pmax ~ 1 leads to an abrupt in¬ 
crease in population diversity from ~ 1 to ~ N. The diversity 
subsequently decreases oc exp (—t/N) (dash-dotted line) 



FIG. 3. Time-aggregated population size distribution. 

The grey shaded area shows the time-aggregated distribution 
of population sizes in our model with 7 = 10“ 9 and N = 1000 
collected over 20 millions collapses. The green and red lines 
show the population size distributions collected, respectively, 
at the very end of each wave and at the very beginning of the 
next wave correspondingly as described in the text. Note that 
they roughly correspond to two peaks of the time-aggregated 
distribution. 

populations immediately after the dominant species with 
Pmax > 1 — l/N was eliminated, but before the available 
biomass was redistributed among all species. At those 
special moments, happening only once every t wave time 
steps, most population sizes are between 7 and 7 -N while 
a small fraction reaches all the way up to ~ l/N. During 
the rapid growth phase immediately after our snapshot 
was taken, all populations grow by the same factor 1/(1 — 
Pmax ) — tV thereby moving all of them to the upper peak 
of the time-aggregated distribution thereby starting the 
new wave. The red curve corresponds to the snapshot 
of all populations immediately after this rescaling took 
place. It shows that at the very beginning of the new 
wave local populations are broadly distributed between 
~ N • 7 and 1 with a peak around l/N. 

Figure 0K shows time-aggregated distributions of pop¬ 
ulation sizes for 7 = 10“ 10 and different values of N rang¬ 
ing between 100 ad 10, 000 while Figure 0)3 shows time- 
aggregated distributions with N = 1000 and for a wide 
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FIG. 4. Time-aggregated distributions for different 
values of N and 7 . Time-aggregated distributions of popu¬ 
lation sizes in our model with A) 7 = 10 -10 and N = 100 
(blue), N = 1000 (red), and N = 10,000 (green). B) 
N — 1000 and varying 7 ranging between 10 -4 (green) to 
10 “ 10 (red) in ten-fold decrements. Note the emergence of a 
nearly universal scale-free tail of the distribution fitted with 
t — 1.7 (dashed line). 


range of 7 . One can see that for 7 < 0 . 01 /AT, the tail of 
the distribution for most abundant populations between 
1/N and 1 is well fitted by a power law 7 r(P) oc 1/ P r_1 ~ 
l/p °.7 (daghgd ii ne i n Figure |4£>) corresponding to the 
power law distribution of population sizes on the linear 
scale dPiob(Pi(t) > P)/dP ~ 1 /P T ~ 1/P 1 - 7 . Overall 
Figures [4]A,B demonstrate that the exponent r for differ¬ 
ent values of 7 and N is remarkably universal. Indeed, 
for a range of parameters where the upper and the lower 
peaks of 7 r(P) are clearly separated, r approaches a uni¬ 
versal value r = 1.7. 

An insight into the origins of the scale-free tail of the 
distribution of population sizes is gained by considering 
a simplified version of our model in which at the start 
of each wave the populations of all species are artificially 
set to be equal to each other resulting in the maximal di¬ 
versity. We further assume that 7 <C 1/N. The passage 
of time t elapsed since the beginning of the current wave, 
leads to a decrease in the number of surviving species 
N SU rv(t ) = Afexp(— t/N), which all have the same pop¬ 
ulation size P = 1/N surv (t) jointly filling up the car¬ 
rying capacity of the environment. Above we ignore a 
negligible fraction (^ 7 ) of the total population of the 
lower peak species. The time-aggregated probability for 
a species to have a population size P{> P — 1/N surv (t) 


is naturally given by N surv (t)/N oc 1/P and thus 

Prob (Pi > P) oc 

o 1 ,/D D 1 dProb(Pj > P) 1 

Prob (Pi = P) =- — -oc 

The exponent r = 2 predicted by this simplified model 
is realized in our actual model for moderately high 7 ^ 
0 . 1 , whereas smaller values of 7 give rise to a different 
universal exponent r ~ 1.7. The decrease of the exponent 
r from 2 to 1.7 in our original model is the result of 
unequal population sizes at the beginning of a new wave. 
In fact, we verified numerically that r = 2 is recovered 
if at the start of each wave one equilibrates all species 
abundances to 1/N. The first section of the SI Text in 
supplementary materials provides additional details on 
how the reduced population diversity D = 1/ Y/ Pf < N 
at the start of population waves affects the exponent r. 

Two panels in Figure [5] illustrate the difference be¬ 
tween the simplified (panel A) and the real (panel B) 
models. In both versions of the model the average jump 
in the logarithm of surviving populations grows exponen¬ 
tially with time t elapsed since the start of the current 
wave: - log(l - Pcollapsed^)) ~ exp(t/N)/N. However, 
unlike the simplified model, the population distribution 
in our real model has a rich hierarchical structure with 
multiple sub-peaks in some waves (color bands in Fig¬ 
ure EP). Remarkably this multi-modal distribution reap¬ 
pears in subsequent waves, implying that the memory 
about the hierarchical structure in the upper part of the 
distribution is transmitted to emerging populations in 
the lower part with sizes starting at 7 . At the start of 
the next wave these same populations will move to the 
upper part of the distribution thereby transmitting the 
history across waves. Colors of symbols in Figure [5] illus¬ 
trate the origin of multiple peaks. Indeed, populations 
in each of these peaks were born during the previous 
wave under similar conditions (the number of substan¬ 
tial populations) as described in the caption. Thus, the 
broadening of time-aggregated population distribution in 
our model compared to its history-free version is a simple 
manifestation of a complex interplay between ” upstairs” 
and ” downstairs” subpopulations transmitting memory 
across waves. 

The population distributions in both upper and lower 
peaks in our model are described by the same exponent 
r. This similarity reflects the fact that individual popu¬ 
lations in the lower peak are exposed to the same multi¬ 
plicative growth as the ones in the upper peak. Finally, 
the intermediate region of the distribution connecting 
two peaks is shaped by rescaling of all populations in 
the lower peak as they are moved up at the beginning 
of a new diversity wave. When peaks are well separated 
(as e.g. for low values of 7 in Figure 0]) the slope of the 
logarithmic distribution in this region has almost exactly 
the same value r — 1 = 0.7 and the opposite sign to the 
slopes in both the upper and the lower peaks. 
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FIG. 5. Memory of population size distribution is pre¬ 
served across several diversity waves. Time course 
of jumps - log[l - Pcoiiapsedit)] in the logarithm of surviving 
populations following a collapse of a substantial population 
Pcoiiapsedit) > 10“ 10 in A) the simplified model in which at 
the start of each wave all populations are set equal to each 
other; B) our basic model. Both were simulated at N = 1000 
and 7 = 10 -20 . Note that our basic model, unlike its sim¬ 
plified counterpart, preserves memory of population sizes dis¬ 
tribution across several subsequent diversity waves. This is 
manifested in similar fractal structure of jumps sizes in waves 
#2-6 shown in panel B). Colors of symbols represent the log 10 
of the number of substantial populations during the the pre¬ 
vious wave, when a given population originated at the small 
size 7 . Thus red dots mark populations originated at the very 
end of the previous wave, while yellow dots - those originated 
when there were two large populations left in the previous 
wave. Finally, green, blue, and purple dots refer to older pop¬ 
ulations in the previous wave. 


Discussion 

In this paper we explore the population dynam¬ 
ics in saturated environments driven exclusively by 
near-complete collapses of sub-populations of competing 
species. This type of dynamics strongly contrasts the 
gradual changes implied in for example the “community 
drift” neutral models [E[ in ecology, or for the most part 
incremental random walks of stock values of individual 
companies. Conversely, collapse-driven dynamics repre¬ 
sents a sudden and usually large change of system com¬ 
position. In ecology such collapses may be caused e.g. 
by invading predators or diseases , whereas in the econ¬ 
omy, companies of any size routinely go bankrupt e.g. 
through excessive debt amplifying the effects of external 
perturbations. 

First, let us consider biological systems. One of the 
predictions of our model is a multimodal logarithmic dis¬ 
tribution of population sizes. Indeed, while the time- 
aggregated distribution is bimodal with distinct upper 


and lower peaks, populations within any given diversity 
wave cluster together in several smaller peaks persisting 
over several waves (color stripes in FigureEB)- This over¬ 
all finding is supported by a growing body of literature 
UMl where multi-modal Species Abundance Distribu¬ 
tions (SAD) in real ecosystems were reported for plants, 
birds, arthropods 0 , marine organisms including single 
cells, corals [12j|, nematodes, fishes, entire seafloor com¬ 
munities 0 , and even extinct brachiopods 0. Like in 
our model, the empirical SADs range over many orders 
of magnitude with a noticeable depletion (or several de¬ 
pletions) at intermediate scales. The magnitude of this 
dip is usually somewhat less than predicted by our ba¬ 
sic model but is consistent with several of its variants 
described below. This includes the model variant #1 in¬ 
spired by the neutral theory of biodiversity [ 5 ] thought 
to apply to a variety of ecosystems including microbial 
communities [6, 7] (see SI Figure in supplementary ma¬ 
terials). 

Needless to say, our model is not unique in generat¬ 
ing multimodal distributions (see e.g. jl3j for other ex¬ 
amples). Conversely, some of the variants of our model 
give rise to interesting population dynamics including di¬ 
versity waves even without any depletion in the middle 
of the log-binned SAD. We argue that a more reliable 
characterization of underlying dynamical processes can 
be obtained from time-series data. First, all systems 
capable of diversity waves are described by rapid large 
changes in populations of individual species. Such sud¬ 
den, population-scale shifts can occur e.g. due to intro¬ 
duced diseases or the disappearance of keystone species 
00 thereby changing the composition of the entire 
food-web. On the microbial scale, sudden invasion of a 
new bacteriophage may lead to multiple orders of magni¬ 
tude reduction in the population of susceptible bacteria 
@, LL8] 2 potentially leading to their complete local extinc¬ 
tion 9[. Phage-driven collapses do not spare bacteria 
with large populations and may even be biased towards 
such as postulated in the Kill-the-Winner (KtW) hypoth¬ 
esis [19]. The magnitude and characteristic timescales of 
population changes in microbial ecosystems is still be¬ 
ing actively discussed in the literature. While Ref. [20] 
reports that over half of all bacterial species in marine 
environments varied between abundant and rare over a 
three-year period, other studies 0 found more modest 
variability at the level of species or genera over weeks 
to months period. However, everyone seem to agree on 
dramatic and rapid (often on the scale of days [22j) pop¬ 
ulation shifts at the level of individual bacterial strains 
@, |[l, [23] caused by phage predation [22]. Except for 
interchangeable gene cassettes (metagenomic islands) re¬ 
sponsible for either phage recognition cites or alterna¬ 
tively resistance to phages [24] these strains routinely 
have very similar genomes and thus may have near iden¬ 
tical growth rates. Hence, they are capable of coexistence 
in the saturated state implicitly assumed in our model. 

Extinctions and collapses in our model are assumed to 
be caused exclusively by exogenous effects such as natural 
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catastrophes or predation by external species not shar¬ 
ing the carrying capacity of our environment. Real-life 
ecosystems can also collapse due to endogenous effects, 

i. e. internal interactions between species. Such intrinsic 
collapse mechanisms were the sole focus of earlier models 
by us and others (see e.g. |9, 25, 26]). 

On vastly longer, geological timescales, the collapse- 
driven dynamics of our model resembles that of species 
extinctions and subsequent radiations in the paleonto¬ 
logical record nm. One example is the recolonization 
by mammals of a number of ecological niches vacated 
by dinosaurs after the end-Cretaceous mass extinction 
thought to be preceded by a gradual depletion of diver¬ 
sity among dinosaurs who were finally wiped out by a 
singular catastrophic event [29]. Interestingly, the ex¬ 
tinction rate of taxonomic orders appears to be indepen¬ 
dent of their size quantified by the number of genera they 
contains 0 , which is also one of the assumptions of our 
basic model of collapse-driven dynamics. 

The second area of applications of our model is to de¬ 
scribe company dynamics in economics. The size or the 
market share of a publicly traded company reflected in its 
stock price is well approximated by a random walk with 
(usually) small incremental steps [30|. However, as in the 
case of ecosystems, this smooth and gradual drift does 
not account for dramatic rapid changes such as bankrupt¬ 
cies or market crashes. In case of companies the main 
driver of sudden changes is their debt [3l| . When the in¬ 
trinsic value of a company is much smaller than its debt, 
even small changes in its economical environment can 
make it insolvent not sparing even the biggest companies 
from bankruptcies (think of Enron and Lehman Broth¬ 
ers). Empirically, the year-to-year volatility of comp any ’s 
market share varies with its size S', A S/S oc S -0,2 [32| . 

Abundance distributions in our original model and 
many of its variants are characterized by a power-law tail 
with an exponent r close to 2. This is in approximate 
agreement with the empirical data on abundance distri¬ 
butions of bacteria in soil samples [33|, marine viruses 
(phages) [3^ . 

In the economics literature, the distributions of com¬ 
pany sizes [35|, as well as those of wealth of individuals 
|36j are known to have similar scale-free tails. Recent 
data for companies [35j and personal wealth [36j sug¬ 
gest 1/P 1,8 tail of the former distribution and 1/P 2,3 
tails of the latter one. Traditionally, scale-free tails in 
these distributions were explained by either stochastic 
multiplicative processes pushed down against the lower 
wall (the minimal population or company size, or wel¬ 
fare support for low income individuals) [37H39| , by vari¬ 
ants of rich-get-richer dynamics [40] , or in terms of Self- 
Organized Criticality [25|, |4l|. The emphasis of the latter 
type of models on large system-wide events (avalanches 
[ 251 .0 or collapses [42| ) and on separation of timescales 
is similar in spirit to the collapse-driven dynamics in our 
models. A potentially important socio-economic impli¬ 
cation of our model is that during each wave contingent 
sub-peaks in the “upstairs” part of the distribution are 


imprinted on the “downstairs” part and thereby can be 
repeated in the new wave following the “revolution”. 

Needless to say, our models were simplified in order 
to compare and contrast the collapse-driven dynamics to 
other mathematical descriptions of competition in satu¬ 
rated environments. The SI Text in supplementary ma¬ 
terials describes several variants of our basic model that 
in addition to population collapses include the following 
elements: 

1. “Neutral drift model” assumes changes of popula¬ 
tion sizes during time intervals between collapses 
as described in Ref. [E[. In this model in addition 
to collapses a population of size Pj r andomly drifts 
up and down A Pi oc ±^/Pi( 1 — Pi). The resulting 
diversity waves and time-aggregated distributions 
can be found in the supplementary SI Figure. 

2 . “Exponential fluctuations model” is another vari¬ 
ant of the neutral scenario where the population 
sizes between collapses undergo slow multiplicative 
adjustments A Pi oc ±QiPi restricted by the overall 
carrying capacity of the environment. Details and 
the resulting time-aggregated distribution can be 
found in the supplementary S 2 Figure. 

3. “Interconnected environments model” is another 
neutral variant of our basic rules in which spa¬ 
tially separated sub-populations of the same species 
are competing with each other for the same nu¬ 
trient. Sub-populations are connected by the dif¬ 
fusion, that is much slower than the diffusion of 
shared nutrient. In this model a collapsed sub¬ 
population is replenished by a small number 7 of 
individuals diffusing from other environments, see 
the supplementary S3 Figure. 

4. ”Kill-the-Winner (KtW) model” where collapse 
probability c systematically increase with the pop¬ 
ulation size as suggested by the studies of phage- 
bacteria ecosystems [l9]. In this particular case 
the diversity dynamics and the scale-free tail of the 
population distribution becomes sensitive to the ex¬ 
tent that the large populations are disfavoured by 
collapse. When the collapse probability is propor¬ 
tional to population size, one obtain a flat distri¬ 
bution where numbers of species in each decade of 
population sizes are equal to each other, see the 
supplementary S4 Figure. 

5. ”Kill-the-Looser (KtL) model”, where collapse 
probability c systematically decreases with the pop¬ 
ulation size P as c(P ) ^ P ~ 0,2 as suggested by 
the empirical studies of company dynamics |32j . 
As seen in the supplementary S5 Figure the diver¬ 
sity dynamics and the scale-free tail of the popula¬ 
tion distribution are both remarkably robust with 
respect to introduction of size-dependent collapse 
rate. 
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6 . ”Fitness model” in which each of the species has 
its own growth rate f^ during rapid re-population 
phase and its own collapse probability <7 . The sup¬ 
plementary S 6 Figure show that the overall shape of 
the time-aggregated distribution is similar to that 
in our basic model, whereas its lower panel illus¬ 
trate the interplay between population size and the 
the two variables that define the species’ fitness. 

7. “Resilience model” as a variant of the above fitness 
scenario, in which collapsing species do not neces¬ 
sarily go into extinction. Instead, each species is 
assigned its own “survivor ratio” 7 $ defined by the 
population drop following a collapse: Pi —» 7 ^. 
As in the previous variant each of the species is 
also characterized by its own growth rate The 
supplementary S7 Figure shows that for interme¬ 
diate populations the time-aggregated distribution 
is described by a power law scaling. Compared to 
the basic model it has a broader scaling regime and 
larger likelihood to have most of the “biomass” col¬ 
lected in one species. 

Captions to supplementary S1-S7 Figures provide more 
detailed description of model dynamics in each of these 
cases. Overall, the ‘ simulations of the variants of our 
basic model described above preserve the general pat¬ 
terns of collapse-driven dynamics such as diversity wave 
dynamics, and a broad time-aggregated distribution of 
population sizes with scale-free tail for the most abun¬ 
dant species. 
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FIG. 6 . Average population vs species’ properties in 
the “Fitness model” variant # 6 . Time-averaged popu¬ 
lation of a species (see color scale on the right) plotted as a 
function of its re-population growth rate Qi (x-axis) and pop¬ 
ulation drop after collapses 7 i (y-axis). in a variant of our 
model with fitness differences between species. Note that the 
population increase with both IT and 7 *. Populations and 
fitness parameters of N — 1000 species were taken from 50 
million snapshots of the model. 

The classic definition of the fitness of a species in terms 


of its long-term exponential growth rate 43] is clearly in¬ 
appropriate for our model. Indeed, the long-term growth 
rate of each of our species is exactly zero. One must keep 
in mind though that fitness is a very flexible term and 
has been defined in in a variety of ways HU] reflecting 
(among other things) different timescales of growth and 
evolution [25]. and relative emphasis on population dy¬ 
namics vs. risk minimization (45j- An appropriate way 
to quantify species’ success in a steady state system like 
ours is in terms of their time-averaged population size 

m)) t . 

In the last two variants of our basic model we add 
fitness-related parameters to each of the species: its 
short-term exponential growth rate f^ (model 6 and 7), 
its propensity to large population collapses quantified by 
their frequency <7 (model 6 ), and the severity of col¬ 
lapses quantified by surviving fraction 7 \ of the popula¬ 
tion (model 7). Fig. 6 shows the average population size 
(. Pi(t))t as function of fb and 7 $ in the model 6 . As one 
can naively expects species with larger short-term growth 
rates and larger surviving ratios also tend to have sub¬ 
stantially larger populations (the red area in the upper 
right corner of Fig. 6 ). 

While in our models the probability <7 and severity 
1 / 7 i of collapses are assumed to be independent of growth 
rate Qi in reality they are often oppositely correlated. 
Indeed, in biology much as in economics/finance fast 
growth usually comes at the cost of fragility and exposure 
to downturns forcing species to optimize trade-offs. 

Some environmental conditions favor the fast growth 
even at the cost of a higher risk of collapse, whereas other 
could call for sacrificing growth to minimize probability 
or severity of collapses. Species in frequently collapsing 
environments considered in our study are known to em¬ 
ploy bet-hedging strategies [ 2 |, l45l447| to optimize their 
long-term growth rate. This is obtained by splitting their 
populations into “growth-loving” and “risk-averse” phe¬ 
notypes [45l, H3, Eg. One example of this type of bet¬ 
hedging is provided by persistor sub-populations of some 
bacterial species consisting of 7 ^ 10 _4 of the overall 
population [49|, [50| increasing to as much as 7 = 10 ~ 2 
in saturated conditions (S. Semsey, private communi¬ 
cations). A bet-hedging strategy with persistor sub¬ 
population 10 _2 somewhat reduces the overall growth 
rate (only by 1 %) while dramatically reducing the sever¬ 
ity of collapses caused by antibiotics. From Fig. 6 one 
infers that this is indeed a good trade-off. 

In this study we presented a general modeling frame¬ 
work for systems driven by redistribution of rate-limiting 
resources freed up by episodic catastrophes. In spite 
of their simplicity the population dynamics in such sys¬ 
tems happens on at least four hierarchical timescales. 
At the shortest timescale the populations grow expo¬ 
nentially repopulating resources vacated during a catas¬ 
trophic extinction event. This exponential growth re¬ 
sults in a steady state at which the system is poised ex¬ 
actly at the carrying capacity of the environment. At 
even longer timescales the system is described in terms 
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of diversity waves that are the main focus of this study. 
These waves are an emergent dynamical property of the 
system in which population of surviving species grow ex¬ 
ponentially while species diversity decays. Remarkably 
the information about the “upstairs” and “downstairs” 
population peaks survives the “revolution” at the end of 
each wave. This memory gives rise to the final, longest 
timescale in our system correlating several consecutive 


waves. All of this complexity is already present in our 
basic one-parameter model. In spite of its simplicity this 
model and its variants provide the foundation for future 
studies of collapse-driven dynamics in ecosystems, mar¬ 
ket economies, and social structures. 
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SUPPLEMENTARY MATERIALS 
Fokker-Planck equation for the basic model 

Let’s consider a version of our model with a very low 
value of 7 to ensure the complete separation between the 
upper and lower peaks. The multiplicative dynamics of 
surviving populations in the upper peak is described by 
the following elementary step: 


m 


m -1) 

1 — Pj(t) (t) 


. Here Pj(t)(t) is the population at time step t of the 
species j(t) that went extinct at this time step. It can be 
also easily integrated for all times since the beginning of 
the current wave at time step 1 : 


m = 


m) 

1 - Et'=i- P i(t')( 1 ) 


(SI) 


Indeed, Et'=i P j(t') ( 1 ) is the total initial (at time step 
1 ) populations of all species that went extinct by the 
time step t. Hence 1 — Ylt'=i ^j 0')(1) " the total 
initial population of all surviving species used to nor¬ 
malize their initial populations to give their popula¬ 
tions at the time step t (population ratios of surviving 
species are preserved in our basic model). Taking the 
logarithm of both sides of Eq. [Si] and approximating 

- iog(i - Eki p i(t')( 1 )) - Eki - p i(t')( 1 )> which holds 

as long as the system is still far away from the end of the 
wave (Et'=i p j(t')(l) ^ r one S ets: 


log Pi(t) = iog p i(i) + Y p m( l ) • ( S2 ) 

t'=1 

The stochastic dynamics within a single wave can thus 
be described by the following equation: 

d\og(Pi{t)) _ 

— j t — - p my l > • ( b3 ) 

The exponent r of the population size distribution in 
our model is determined by the balance between the noisy 
multiplicative population dynamics and the exponential 
loss of surviving species due to collapses. It can be ap¬ 
proximated by the following Fokker-Plank-like equation: 

dn dn 8 2 tt 

m = ^ p m^ + *'^^ +eplsodlc source terms 

(S4) 

Here 7r(logP, t) is the time-dependent population abun¬ 
dance distribution in the upper peak, p - the loss term due 
to population collapses, p - the logarithmic drift velocity 
and a - is the logarithmic dispersion, which is totally ab¬ 
sent in the simplified model where all populations start 
at the same size. 

The Eas lS2llS3l allow us to derive the parameters of 
the Fokker-Plank equation in terms of the distribution of 


population sizes P^(l) at the start of the wave. Indeed, 
early in the wave one has p = 1/7V, p = (Pj(l))^ = 1/N 
and o’ = W l) 2 ), - <P;(1 )} 2 = (1/N) ■ (1/D(1) - 1/N). 
Note an unusual connection between the population di¬ 
versity at the start of a wave D( 1), and the diffusion 
coefficient a in the Fokker-Plank equation. 

The stationary solution for the time-aggregate distri¬ 
bution f 7 r(log P, t)dt has an exponential tail exp( —(r — 
l)logP). It corresponds to the power law tail of the 
species population distribution oc P -r . The exponent r 
is defined by one of the two solutions to the quadratic 
equation 


0 = --I + -I(T-l) + cr(T-l) 2 . (S5) 

In the version of the model, where all populations at the 
start of the wave are equal to each other, sizes of surviv¬ 
ing populations increase deterministically as exp (t/N)/N 
(see main text for the derivation) and thus have zero dis¬ 
persion: <j = 0. Hence in this simplified version the 
exponent r = l + p//i = 2is determined by balancing 
only the first two terms of this equation. 

We have numerically verified that the decrease of the 
exponent from r = 2 in the simplified model down to 
r = 1.7 in our original model is driven entirely by noise 
(unequal population sizes) resulting in a finite value of 
a. A non-zero value of a in the Eq. [S5] results in r < 2. 
For example, if the populations at the start of each wave 
had a Poisson distribution so that a = p = 1/7V, the 
exponent r = (\ + y/h)/2 ce 1.62 would have been defined 
by the solution of the golden mean equation 0 = —1 + (t — 
1) + (r — l) 2 . While currently we have no first-principles 
argument allowing us to derive the value of a in our basic 
model, the result from a Poisson distribution is not too 
far from the empirically observed exponent r = 1.7 


Model variants 

To test the robustness of our basic model with respect 
to rule changes we considered the following seven vari¬ 
ants: 

1. “Neutral drift model”. 

This variant extends our basic model by adding 
to our standard model the random neutral drift 
of population sizes (Hubbell SP (2001) The uni¬ 
fied neutral theory of biodiversity and biogeog¬ 
raphy (MPB-32), Princeton University Press) be¬ 
tween subsequent collapse events. To simulate 
this random drift, at every time step the popula¬ 
tion of each species changes up or down as pre¬ 
scribed by dis persion of bi nominal distribution: 
Pi —>> Pi dz yjr • Pi(l — Pi ), where r is the pa¬ 
rameter quantifying the magnitude of fluctuations 
proportional to the inverse of the total population 
size and the square of the birth/death rate. Af¬ 
ter drift changes were applied to all populations 
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we rescale them back to their carrying capacity 
^2 Pi = 1. This is followed by a collapse event as 
in the standard model. Fig. [SU illustrates typical 
time courses of the diversity D(t) = 1 /^2Pi{t) 2 
and time-aggregated species abundance distribu¬ 
tions in this model variant for three values of r and 
compares them to our basic model. 



20000 30000 

time (collapse events) 



Iog10(population size) 

FIG. SI. “Neutral drift model”. This variant extends our 
basic model with N = 1000 and collapse ratio 7 = 10 -9 by 
adding the neutral drift at rate r taking place between sub- 
se quent collapse events in our standard model: Pi Pi d= 
yfr - Pi(l — Pi). The lower panel shows the time-aggregated 
distributions in our system simulated for 10 6 collapse events. 
The grey shaded area refers to our basic, unmodified model, 
i.e. to the r — 0 case, while three color lines correspond to 
r — 10 -9 (blue), r — 10 -7 (green), and r — 10 -5 (red). 
The upper four panels illustrate typical time courses of the di¬ 
versity D{t) = 1 /^2Pi(t) 2 in our basic model and for three 
values of the r color-coded as in the lower panel. 


2. “Exponential fluctuations model”, where pop¬ 
ulations are exposed to random exponential shifts 
between successive collapse events. In this ver¬ 
sion of the model between successive collapse events 
all populations exponentially shift thus adjusting 
the way carrying capacity is divided between the. 
This model is similar to the ” Neutral random drift” 
model ffl above except that changes are propor¬ 
tional to Pi and not to \fPi. The population of 
each species is characterized by its own exponen¬ 
tial rate Gi(t) given by n • randi(t ), where n quan¬ 
tifies the overall rate of the redistribution, while 
randi (t) - a random number uniformly distributed 
between 0 and 1 - represents species- and time- spe¬ 
cific shifts. We assume that differences in growth 
rates between species are not constant but instead 
fluctuate on the timescale when a single species col¬ 
lapses. Hence, after each collapse event we reset 
the exponential growth we randomly reset the rates 
Gi(t ) for all species (and not just of the collapsed 
one). In this version of the model we also take into 
account the stochastic nature of the time interval 
r r (t) between two successive collapse events, which 
is randomly chosen from the exponential distribu¬ 
tion with mean value equal to 1. Thus between two 
successive collapse events each species population 
changes as Pi Pi • e G *d) r r(t) anc [ subsequently 
rescaled to the carrying capacity of the environ¬ 
ment: ^2 Pi — 1. As in our basic model, at every 
time step one randomly selected population i col¬ 
lapses and is reset to Pi 7 while all other popu¬ 
lations are rescaled to fill up the carrying capacity: 
"22 Pi = 1. Fig. IS2l shows the simulations of this 
model for several value of n compared to our basic 
(n — 0 ) model. 
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FIG. S2. “Exponential fluctuations model” with N = 
1000 , 7 = 10 -9 , and n = 0.02 (green), n — 0.1 (blue), n — 
0.5 (red) system simulated for 10 6 collapse events. The grey 
shaded area shows the time-aggregated population distribution 
in our basic model, corresponding to the n — 0 limit. 


3. “Interconnected environments model” with 
diffusion. In this version of the model there is a 
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single species distributed between N local envi¬ 
ronments. As in our basic model at every time 
step one local population i is selected for collapse 
and reset to Pi —>• 0 after which the populations 
are normalized back to the carrying capacity. The 
diffusion takes a fraction 7 of the total population 
of 1 and distributes it equally between all local en¬ 
vironments: Pi —>• Pi(l — 7 ) + j/N. Note, that 
here we implicitly assume that populations in all of 
these environments share the same carrying capac¬ 
ity. This is the case when diffusion rate of the rate- 
limiting nutrient is much faster then that of popula¬ 
tions themselves. The main difference of this model 
from earlier variants is that populations in the lower 
peak with Pi <C 1/N grow approximately linearly 
in time (as opposed to exponentially in other ver¬ 
sions of the model). The rate of this linear growth 
is the same for all species and is equal to 'y/N. The 
exponential growth is restored for populations that 
are larger than average. Fig. S3 shows the time- 
aggregated distribution of local population sizes in 
this model variant. 



Iog10(population size) 


ulation size. At each time step we select a ran¬ 
dom population to collapse with probability oc P?. 
As before the collapsing population is reset to 
Pi ^ 7 and all populations are subsequently grown 
with equal exponential rates to complete satura¬ 
tion: ^2 Pi — 1. Fig. S4 examines time-aggregated 
population distributions in KtV model variant for 
different values of a. Whereas small and moder¬ 
ate a preserve diversity wave dynamics, the a = 1 
version does not exhibit diversity waves and pre¬ 
dicts a population size distribution distributions, 
dP/ds oc 1 /s, or dP/d\og(s) = constant (equal 
number of species in each decade of population 
sizes). 



Iog10(population size) 

FIG. S4. “Kill-the-Winner (KtW) model” in which 
larger populations are preferentially targeted for collapse: a oc 
Pi . Different colors correspond to time-aggregated SADs in 
the model with N — 1000, 7 = 10 -9 , and a = 0.01 (green), 
0.2 (blue), and a = 1.0 (red) simulated over 5 • 10 6 collapse 
events. The grey shaded area refers to time-aggregated pop¬ 
ulation distribution in our basic, unmodified model with the 
same N and 7 . 


FIG. S3. “ Interconnected environments model”. Fig¬ 
ure show time-aggregated species abundance distributions in 
the model with N — 1000 environments connected by diffu¬ 
sion of strength 7 = 10 -9 simulated over 10 6 collapse events. 
The standard model with the same parameters is shown as the 
grey shaded area. 

4. “Kill-the-Winner (KtW) model”. For bacte¬ 
rial populations the direction of the trend (if any) 
of collapse probability with population size is cur¬ 
rently unknown. In fact one can plausible make 
a case for increasing of the probability of collapse 
with population size due to larger populations mak¬ 
ing easier to find and overall and more attractive 
targets for phages. In microbiology preferential 
targeting of large bacterial populations by viru¬ 
lent phages is known as ” Kill-the-Winner” (KtW) 
hypothesis (Thingstad TF and Lignell R (1997) 
Aquatic Microbial Ecology 13:19-27). Here we sim¬ 
ulate the version of our basic model where the col¬ 
lapse probability systematically increases with pop- 


5. “Kill-the-looser (KtL) model” in which a col¬ 
lapse probability declines with population size in 
a power-law fashion. At each time step one select 
a random population to collapse with probability 
oc P ~°‘ 2 where Pi is the current population size of 
species i. In economics this corresponds to an in¬ 
tuitively plausible notion that larger companies are 
less likely to go bankrupt than smaller ones. Em¬ 
pirically, this trend is described by a power law with 
the exponent - 0.2 (Nunes Amaral LA, Buldyrev 
SV, Havlin S, Leschhorn H, Maass P, Salinger MA, 
et al. (1997) Journal de Physique I. 7: 621-633.) 
Notice the emergence of the lower peak distribu¬ 
tion distribution above 7 that is much more nar¬ 
row than in our standard model. This makes sense 
as in the course of each wave small populations 
tend to collapse over and over. These repeated col¬ 
lapses don’t drive other populations up by much 
and thus their only consequence is clustering of 
small populations close to the very bottom of the 
lower peak distribution at and above 7 . When the 
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dominant upper peak population finally collapses 
all small populations are rescaled up to form a nar¬ 
row distribution around 1/AT. This is very similar 
to our simplified memory-free model described in 
the main text and shown in Fig. 4A. The new 
wave starts with very high diversity D(t) ~ N 
which is subsequently reduced with time as D(t) = 
N sur v{t) — Nexp(—t/N). Here we ignore a rela¬ 
tively small AT 0,2 -fold decline in collapse frequency 
over the range of population sizes between 1/N and 
1. Within the same approximation each surviv¬ 
ing population grows as P(t) oc exp (t/N). The 
time-averaged distribution of populations thereby 
approaches the scaling regime described by: 


Prob (Pi(t) > P ) 


1 

p ■ 

dProb(Pj (t) > P) 1 
dP ~ P^ 


(S 6 ) 


In reality the scaling exponent of the tail is around 
1.8. It is the same as in our standard model but 
for a different reason. Indeed, taking into account 
that lifetime of a population before collapse scales 
as 1/P - 0 - 2 = P 0 - 2 one gets Prob(P i (t) = P) - 

p °- 2 l 

p 2 pi.8 • 


1000 


~ 100 

co 
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time (million collapses) 
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FIG. S5. “Kill-the-loser (KtL) model” in which the col¬ 
lapse probability declines with population size. The figure 
shows an N = 1000, 7 — 10 -9 system simulated for 10 6 
collapse events. The upper panel illustrates the recurrent di¬ 
versity waves, whereas the lower panel shows time-aggregated 
distributions, with the grey shaded area referring to our stan¬ 
dard model. 


6. “Fitness model” with heterogeneous, species- 
specific growth rates and extinction probabilities. 


Each species is assigned a growth rate used when 
it repopulates the freed-up carrying capacity of the 
environment. It also has its own extinction prob¬ 
ability C{. Both Tli and <7 are logarithmically dis¬ 
tributed in the interval between 0.1 and 1. That 
is to say their log 10 are uniformly distributed be¬ 
tween — 1 and 0 . At each time step we select one of 
N populations, with probability oc c&, this species 
goes extinct. It is immediately replaced by a new 
species with the population P k j = 10 -9 , new 
growth rate and extinction probability < 7 . Sub¬ 
sequently all of the populations i = 1,2...TV are 
rescaled proportional to their growth rates 


to fill up the carrying capacity of the environment 
^2 Pi = 1. The upper panel in Fig. S6 shows that 
the time-aggregated population distribution in this 
model preserves its power-law tail, whereas lower 
panel illustrates that in order for a species to reach 
substantial population size its fitness parameters 
need to be particularly favorable. Indeed, popula¬ 
tions larger than 1/N = 0.001 tend to have smaller 
than average extinction probabilities < 7 , and larger 
than average growth rates 

7. “Resilience model” where heterogeneous, 
species-specific growth rates and survival ratios 
(population reduction following a collapse) are 
competing with each other. Each species is 
assigned a growth rate G [0.1,1] and collapse 
size 7 i G [10 - 9 ,10 -2 ], both logarithmically dis¬ 
tributed (uniform distribution of the logarithm 
of the variable). At each time step we select one 
of N populations, and collapse its population 
Pk lk m Pk- Note that unlike in previous versions 
we scale down the population proportional to its 
size and not proportional to the carrying capacity 
of the environment. This reflects a different 
interpretation away from our basic model, where a 
collapse represents not an extinction of the species 
followed by the appearance of a new species at 
a fixed (very small) population size. In the new 
version of the model a collapse represents a sudden 
but proportionate reduction of a population e.g. 
due to species’ phenotypic or genotypic bet¬ 
hedging. In this version of the model an extinction 
happens only if a very low population is reached, 
i.e. when 7 k Pk < 10 -9 . If this lower bound is 
reached the old species goes extinct and a new 
species with the initial population P k = 10 -9 is 
introduced. The new species is assigned new ran¬ 
dom values of Sl k and 7 ^. As in the previous model 
following each collapse all populations i = 1 , 2 ... AT 
are rescaled proportional to their growth rates 
Pi ^ Pi + (P k - max( 7 fe P fc , 10 -9 )) • G’ to 
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FIG. S6. “Fitness-model” with heterogeneous, species- 
specific growth rates and extinction probabilities. Each species 
is assigned a growth rate Qi used when it repopulates the freed- 
up carrying capacity of the environment. It also has its own 
collapse probability Ci. Both Qi and Ci are logarithmically dis¬ 
tributed in the interval between 0.1 and 1. The purple curve 
in the upper panel shows the time-aggregated population dis¬ 
tribution whereas the grey shaded area is that for the standard 
model where species’ growth and collapse rates are all equal to 
each other. The lower panel shows the average collapse proba¬ 
bility and growth rate binned by the size of the population, he 
average growth rate (fh) (blue) and the average collapse proba¬ 
bility (a) (red shaded area) of species binned by their collected 
at every time step. Both curves represent time-aggregated av¬ 
erages as individual populations change with time. 


fill up the carrying capacity of the environment: 

E p i = !• 
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FIG. S7. “Resilience model” with heterogeneous, species- 
specific growth rates and survival ratios following a collapse. 
Each of the N m. 1000 species is assigned a growth rate Qi G 
[0.1,1] and collapse size^i G [10 —9 ,10 _2 ] ; both logarithmically 
distributed. A) The blue curve shows the time-aggregated pop¬ 
ulation distribution, whereas the grey area refers to that in our 
standard model from the main text. B) The average growth 
rate (Qi) (blue) and the average survival ratio multiplied by 
200 to have the same range in the plot ( 7 f) (red shaded area) 
as a binned by the population size collected at every time step. 
Both curves represent time-aggregated averages as individual 
populations change with time. C) The average (arithmetic) 
population size as a function of species’ survivor ratio 7 
D) The average (arithmetic) population size as a function of 
species ’ growth rate Qi. 













